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Starting from the spectrum of the radially symmetric quantum harmonic oscillator in two di- 
mensions, we create a large set of nonlinear solutions. The relevant three principal branches, with 
nr = 0, 1 and 2 radial nodes respectively, are systematically continued as a function of the chemical 
potential and their linear stability is analyzed in detail, in the absence as well as in the presence 
of topological charge m, i.e., vorticity. It is found that for repulsive interatomic interactions only 
the ground state is linearly stable throughout the parameter range examined. Furthermore, this 
is true for topological charges m = or m = 1; solutions with higher topological charge can be 
unstable even in that case. All higher excited states are found to be unstable in a wide parametric 
regime. However, for the focusing/attractive case the ground state with Ur — and m = can 
only be stable for a sufficiently low number of atoms. Once again, excited states are found to be 
generically unstable. For unstable profiles, the dynamical evolution of the corresponding branches 
is also followed to monitor the temporal development of the instability. 



I. INTRODUCTION 



The study of trapped Bose-Einstein condensates 
(BECs) has had a high impact in recent years in a 
number of fields, including atomic, molecular, and opti- 
cal physics, nuclear physics, condensed matter physics, 
chemical physics, applied mathematics, and nonlinear 
dynamics [J, H, H]- From the point of view of the lat- 
ter, the topic of particular interest here is that at the 
mean field level, the inter-particle interaction is repre- 
sented as a classical, but nonlinear self-action [1], leading 
to the now famous Gross-Pitaevskii (GP) equation as a 
celebrated model for BECs in appropriate settings. This 
has resulted in a large volume of studies of nonlinear ex- 
citations, including the prediction of excited states [5[, 
the experimental observation of dark 0, 0, B Q , bright 
dO, EH and gap [l2[ solitons in quasi-one-dimensional 
systems, as well as theoretical and experimental inves- 
tigations of vortices, vortex lattices [3, [13], and ring 
solitons [III, [IB, in quasi- two-dimensional systems. 

Apart from purely nonlinear dynamical techniques, 
such as the perturbation theory for solitons employed 
in Ref. [l5[, other methods based on the corresponding 
linear Schrodinger problem may also be employed for the 
study of excited BECs. In particular, the underlying lin- 
ear system for a harmonic external trappingpotential is 
the quantum harmonic oscillator (QHO) [H, [1, 0], whose 
eigenstates are well-known since Schrodinger 's original 
treatment of the problem in 1926. Beginning with the 
linear equation, solutions can be numerically continued 
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to encompass the presence of the nonlinear representation 
of the inter-particle interaction. Then, the QHO states 
can persist, bifurcate, or even disappear. This path does 
not seem to have been exploited in great detail in the 
literature. In one spatial dimension, it has been used in 
Ref. [1^, [1^, to illustrate the persistence [l8| and dy- 
namical relevance [HI of the nonlinear generalization of 
the QHO eigenstates. In higher dimensions, the work of 
Ref. [20| illustrated the existence of solutions in a radially 
symmetric setting. Further progress has been hampered 
by the additional difficulties in (a) examining the linear 
stability and (b) converting the radial coordinate system 
to a Cartesian one to study evolution dynamics, includ- 
ing nonlinear stability. However, the mathematical tools 
for such an analysis exist, as we discuss in more detail 
below, and have to a considerable extent been used in 
Ref. [13, especially in connection with ring-like struc- 
tures with vorticity. 

In this paper we present a systematic analysis of the 
existence, linear stability, and evolution dynamics of the 
states that exist in the spectrum of harmonically confined 
condensates from the linear limit, and persist in the non- 
linear problem. Our analysis shows that in the case of 
repulsive interactions (defocusing nonlinearity ) , the only 
branch that is linearly stable consists of the ground state 
solution, i.e., a single hump with = radial nodes. All 
higher excited states are linearly unstable and break up 
in ways that we elucidate below, if evolved dynamically 
in our system; typically this last stage involves also the 
loss of radial symmetry to unstable azimuthal perturba- 
tions. Here, we will treat explicitly the cases of = 0, 1 
and 2 radial nodes; we have confirmed similar results for 
higher number of nodes. 

On the other hand, in the case of attractive interactions 
(focusing nonlinearity), the system is subject to collapse 
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in the absence of the potential. In fact, this constitutes 
the critical dimension for the underlying nonlinear prob- 
lem [m beyond which collapse is possible. In this case, 
we observe that only the ground state branch with small 
values of the chemical potential, i.e., a small number of 
atoms, can be marginally stable. Once again, all excited 
states are unstable and the typical scenario here involves 
the manifestation of collapse-type catastrophic instabili- 
ties [2l|, as we elucidate through numerical simulations. 

Our presentation is structured as follows: in section 
II, we provide the theoretical setup and numerical tech- 
niques. In section III, we present and discuss the relevant 
results. In section IV, we summarize our findings and 
present our conclusions. Finally, in the appendix we pro- 
vide relevant details regarding our numerical methods. 



II. SETUP AND NUMERICAL METHODS 

We will use as our theoretical model the well-known 
GP equation in a two-dimensional (2D) setting. As is 
well-known, the "effective" 2D GP equation applies to 
situations where the condensate has a nearly planar, 
i.e., "pancake" shape (see, e.g., Ref. and references 
therein^ We express the equation in harmonic oscillator 
units [23| in the form: 
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where u = u{x^y,t) is the 2D wave function (the sub- 
script denotes partial derivative). The wave function is 
the condensate order parameter, and has a straightfor- 
ward physical interpretation: p = \u\'^ is the local con- 
densate density. The external potential V{r) = A^r^/2 
assumes the typical harmonic form, with A being the ef- 
fective frequency of the parabolic trap; the latter can 
be expressed as A = ujr/^z^ where ujr,z are the con- 
finement frequencies in the radial and axial directions, 
respectively. It is assumed that A <C 1 for the pan- 
cake geometry considered herein; in particular, we choose 
A = 0.1 for our computations. Finally, cr = ±1 is the nor- 
malized coefficient of the nonlinear term, which is fixed 
to —1 for attractive interactions and to +1 for repulsive 
interactions. Accordingly, the squared norm is 
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where TV is the number of atoms and U = cr\U\ is the 
usual nonlinear coefficient in the GP equation in har- 
monic oscillator units, renormalized to two dimensions 
appropriately [2^ . 

In order to numerically identify stationary nonlinear 
solutions of Eq. ([T]), we use the standing wave ansatz 
u{x^y^t) = v{r) exp[i(/it+m^)], where we assume the 
density of the solution is radially symmetric with chemi- 
cal potential /i and topological charge (vorticity) m. This 



results in the steady-state equation: 
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Equation (|3]) exhibits infinite branches of nonlinear 
bound states, each branch stemming from the corre- 
sponding mode of the underlying linear problem. These 
branches are constructed and followed using the method 
of pseudo-arclength continuation [25|, [2^, in order to ob- 
tain subsequent points along a branch, once a point on 
it has been identified. The first point on each branch is 
found using a bound state of the underlying linear equa- 
tion 
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The linear state corresponds to a solution of Eq. (|3|) with 
parameters fi = frequency A, vorticity m, number of 
radial nodes rir and normalization J 27rr dr\v\^ tending 
to zero; in that limit, the nonlinear ity becomes negligi- 
ble and the linear solutions are the well-known Gauss- 
Laguerre modes with E = {2nr + m + 1)A. Using such 
an initial guess, a non-trivial solution is found through 
a fixed-point iterative scheme for a slightly perturbed 
value of the chemical potential. This is done on a grid of 
Chebyshev points suited to the radial problem, following 
the approach of Ref. for the Laplacian part of the 
equation, as is explained in detail in the appendix. The 
pseudo-arclength method, used to trace subsequent solu- 
tions, works via the introduction of a pseudo-arclength 
parameter s and an additional equation F{v^ /i, s) = 
such that F{v,/2^0) = where (v^fl) is a solution of 
Eq. (|3]). We used for F: 
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Lastly, the new expanded system of equations is solved 
using a predictor-corrector method. 

The linear stability of the solutions is analyzed by using 
the following ansatz for the perturbation: 
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where the asterisk stands for the complex conjugation 
and is the polar angle. One solves the resulting 
linearized equations for the perturbation eigenmodes 
{a(r), b{r)} = {aq(r), bq{r)} and eigenvalues A = Ag asso- 
ciated with them. The key observation that allows one to 
carry this task through is that the subspace of the differ- 
ent azimuthal perturbation eig enmodes, i.e., subspaces 
of different decouple [13, leading each eigenmode 
e*^^ to be coupled only with its complex conjugate. This 
in turn allows the examination of the stability of the 2D 
problem in the form of a denumerably infinite set of one- 
dimensional radial eigenvalue problems that are solved 
on the same grid of Chebyshev points as the original ex- 
istence problem of Eq. (|3|). 
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FIG. 1: Linear stability analysis along the first four branches: 
ground state (no radial nodes) and l^*-3^^ excited states 
(rir = 1,2 and 3 radial nodes) depicted in the curves from left 
to right in each panel. The top, middle and bottom panels 
display a times the norm of the solution as a function of the 
chemical potential for m = (no topological charge), m = 1 
(singly charged) and m = 2 (doubly charged) solutions. Each 
point on the branch is depicted by crosses for stable solutions 
and diamonds for unstable solutions. The areas of stability 
for m = 0, 1 appear only for rir — and only for lower powers 
for (7 = — 1 and all powers of the ground state for a = 1. 
In the case of m = 2, the ground state with rir = is only 
linearly stable for stronger nonlinearity/norm NU in the case 
of cr = 1. 



In order to examine the dynamics of the cases found to 
be unstable, to confirm the validity of the cases identified 
as stable, and to test for nonlinear instabilities, the ID 
radial solutions of Eq. (|3]) were taken as an initial con- 
dition for a code which solved Eq. ([1]) with a Chebyshev 
spectral radial-polar method in space and a fourth-order 
integrator in time. The Chebyshev spectral radial-polar 
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FIG. 2: (Color online) Stability eigenvalues for chargeless 
(m = 0) solutions. The left column of plots corresponds to 
the real part of the primary eigenvalue for ^ = 0, 1, 2, 3, 4 for 
the ground state with rir = (top panel) and first and second 
excited states with rir = 1 and 2 (middle and bottom panels, 
respectively) for a = — 1; while the right column of plots cor- 
responds to the case of a = +1. The ground state for a — 
is omitted since the entire branch is stable. 
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FIG. 3: (Color online) Same as Fig.[2]for singly charged (m ■ 
1) solutions. 
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FIG. 4: (Color online) Same as Fig.[2]for doubly charged (m = 
2) solutions. Notice that in this case the rir — branch can 
be unstable even for the repulsive interactions, i.e., defocusing 
nonUnearity, see top right panel. 




-0.002 
t=0 



t=600 



0.002 -0.1 

t=3 t=0 





t=9 




t=167 



t=250 



FIG. 5: Data for the ground state {ur — 0) for m = for 
a = —1 (left panels) and a = +1 (right panels). The top 
panels show the profile of the solution and the middle panels 
show the corresponding eigenvalues for ^ = 0, 50 in Eq. JS]) 
on the complex plane (Ar,Ai) of eigenvalues A = \r + 
see Table H] for the correspondence between ^-values and the 
different symbols used to depict the eigenvalues. The bottom 
panels depict the time evolution of the solution in harmonic 
oscillator units after an initial random perturbation of 10~^. 
Solutions in this figure correspond to /i = — 2 for a = — 1 and 
norm N\U\ — 100 for a — For all other figures we use 
11 = —0.5 when a = —1. 



method is the most natural choice for our setting, as it a) 
avoids the conversion of the radial solution into an inter- 
polated Cartesian grid, and, more importantly, b) avoids 
spurious effects associated with a mismatch between the 
symmetry of the solution and that of the grid. For ex- 
ample, we have observed that the use of a Cartesian grid 
artificially enhances the excitation of modes that have 
a similar symmetry as the grid, and in particular the 
q = 4 mode. The results of all of the above numerical 
techniques, i.e., existence, linear stability, and dynamical 
evolution, are reported in what follows. 



III. RESULTS 

Our steady state and stability results are summarized 
in Fig. [TJ The panel shows the continuation from the 
linear limit of the relevant states. The nonlinear ity leads 
to a decrease of the chemical potential with increasing 
power for the focusing case and to a corresponding in- 
crease for the defocusing case. Notice that some of the 
branches presented here have also appeared in the earlier 
work of Ref. 0]. However, the important new ingredi- 
ent of the present work is the detailed examination, both 
through linear stability analysis and through dynamical 
evolution, of the stability of these solutions. The latter is 
encapsulated straightforwardly in the symbols associated 
with each branch. The plus symbols denote branches 
which are stable, while the diamond symbols correspond 



to unstable cases. We observe that the only truly stable 
solution corresponds to the ground state with = 
of the defocusing problem, which for large number of 
atoms can be well approximated by the Thomas-Fermi 
state ^j, i2|]. Furthermore, this is true only for topological 
charges m = and m = 1; for higher topological charges 
m > 2, there are regions of stability for = 0, as was 
originally shown in Refs. [l^ and [sO]. All higher excited 
states with one, two or more nodes are directly found 
to be unstable, even close to the corresponding linear 
limit of these states. In the focusing case, the situation 
is even more unstable due to the catastrophic effects of 
self- focusing and wave collapse [U, |3l|. In particular, 
even the ground state with = may be unstable due 
to collapse (represented by the mode with q = in this 
case), although the instability growth rate may be very 
small, as is the case for m = 0, see e.g.. Fig. [2l Excited 
states are always unstable for the focusing nonlinearity 
also; in fact, they are more strongly so than in the de- 
focusing case, again due to the presence of the q = 
mode. 
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TABLE I: Table of symbols used for the unstable eigenvalues 
in Figs. [5HT3I For eigenvalues smaller than 10~^ or ^ > 11 we 
use a small black dot. 

Having offered an outline of the stability properties of 
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FIG. 8: Same as Fig. [5] for m = 1. 



FIG. 11: Same as Fig. [5] for m = 2. 
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FIG. 12: Same as Fig. \TT\ (m — 2) for the first excited state 
{rir = 1). 




FIG. 13: Same as Fig.[TT](m = 2) for the second excited state 

{Ur 2). 



our solutions in Fig. [TJ we now examine the subject in 
detail in the following results. In Figs. [2H31 we depict 
the real part of the primary {q = 0,1,2,3,4) eigenval- 
ues for m = 0, m = 1, and m = 2, respectively and 
for the states with = 0, 1 and 2 in each case, i.e., 
the ground state and first two excited states in the top, 
middle and bottom panel of each figure. The stability 
results showcase the presence of at least one unstable 
eigenvalue. In fact, there is only one instability eigen- 
value for the ground state of the attractive case, and a 
larger number of such eigenvalues for excited states. The 
presence of such eigenvalues, whose larger magnitude is 
tantamount to more rapid dynamical development of the 
instability, indicates that the dynamical evolution of such 
stationary states, when small perturbations are added to 
them, will manifest the presence of such unstable eigen- 
modes through the deformation and likely destruction of 
the initial structure. 

This is shown in detail for m = in Figs. [SHZl In the 
top row of panels in Fig. [5] we depict the radial profiles 



of the solution for /i = — 2 and a = —1 (left panel) and 
N\U\ = 100 and a = +1 (right panel) — for all other 
figures we use /i = —0.5 when a = —1. The middle 
row of panels in the figure depicts their corresponding 
linear stability spectra in the complex plane (Ar,Ai) = 
(Re(A), Im(A)). The unstable eigenvalues for different 
values of q are given different symbols as described in 
the table in Table H Eigenvalues with Re(A) < 10~^ or 
for g > 11 are plotted with a small dark dot. 

Finally, the bottom row of panels in Fig. [5] depicts the 
corresponding time evolution of the chosen profile after 
a random perturbation of amplitude 10~^ was added to 
the steady state profile at time t = 0. For the profile 
considered in Fig.[5l namely the ground state with m = 0, 
it is clear that the solution for the repulsive case (cr = 
+1) is stable since it corresponds to the Thomas- Fermi 
ground state. On the other hand, the solution for the 
attractive case (a = —1) is weakly unstable, due to the 
q = mode, as discussed above. The eigenvalues for this 
mode are depicted by the circles in the middle-left panel. 
This instability is due to the well-known collapse of the 
solution for the attractive case, where the solution is seen 
to tend to a thin spike carrying all the mass (see the time 
evolution in the bottom- left panels). 

In Figs. [6] and [3 we present the equivalent results for 
the first and second excited states (n^ = 1 and 2 respec- 
tively) and for the same m = case. In these cases the 
instability manifests itself with the presence of a richer 
scenario of unstable eigenvalues. Let us explain in detail 
the first excited state in Fig. [6l The first excited state 
in the attractive case (left panels) is still prone to col- 
lapse as evidenced by the quartet of q = eigenvalues 
depicted by the circles in the middle panel. This quar- 
tet is responsible for the collapse of the central spike of 
the solution as seen in the time series evolution (bottom 
panels). Furthermore, the q = 3 and q = 4 modes (^ and 
□ symbols in the figure) are the most unstable ones and 
are responsible for the azimuthal modulational instabil- 
ity of the first ring of the solution as evident in the time 
evolution (bottom panels). 

On the other hand, for the repulsive case (right pan- 
els) of the first excited state with m = it is clear that 
the q = eigenvalue is stable as there is no collapse in 
the repulsive case. The dynamic evolution of the insta- 
bility in this case leads to the competition between the 
unstable eigenvalues ^ = 3, q = 4 and g = 2 (in order of 
strength) that is seen to be dominated by the most unsta- 
ble mode <7 = 3, as can be evidenced by the three humps 
(surrounding the central peak) displayed at t = 150. 

Finally, in Fig. [71 we depict similar results for the sec- 
ond excited state {ur = 2) for m = 0. As can be seen 
from the figures, the more excited the state, the richer 
the (in) stability spectra. It is worth noticing again that 
the attractive case is, as before, prone to collapse due to 
a strong unstable q = mode while, naturally, the repul- 
sive case lacks this q = collapsing mode. The richer set 
of eigenvalues for higher excited states is easy to inter- 
pret since higher excited states possess more radial nodal 
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FIG. 14: (Color online) Three-dimensional iso-density con- 
tour plots. Top plot: iso-density contour for the first excited 
state rir — 1 for cr = +1, m = 2, i.e., the solution corre- 
sponding to the right panel in Fig. [12] Clearly observable is 
the rotation of the unstable q — A mode due to the intrinsic 
vorticity (m = 2) of the solution. Bottom plot: iso-density 
contour for the ground state Ur — for a = — 1, m = 1, 
i.e., the solution corresponding to the left panel in Fig. [H] 
The rotation of the ^ = 3 mode is quickly arrested and the 
three spikes become stationary. The contours correspond to 
surfaces with density equal to half of the maximum density. 



rings. 

In general, for the repulsive cases where collapse is ab- 
sent, we observe that in all cases the solution is subject 
to azimuthal modulations that produce coherent struc- 
tures reminiscent of the "azimuthons" of Ref. [32[. On 
the other hand, in the attractive cases collapse is ubiqui- 
tous, due to the mode with g = 0; in addition, azimuthal 
modulations emerge. 

The above discussion describes in detail the stability 
and dynamics for the non-topologically charged (m = 0) 
solutions. Let us now describe the case when the solu- 
tions have an intrinsic topological charge, namely m > 0. 
In Figs.lSHTOlwe display the results for the ground state, 
first and second excited states (n^ = 0, 1 and 2 radial 
nodes) with unit topological charge (or vorticity), namely 
m = 1. In fact. Fig. [8] corresponds to the singly charged 
vortex solution at the center of the harmonic trap. As it 
is well known, this solution is stable in the repulsive case 
(see right panels) and unstable in the attractive case (see 
left panels). It is interesting to note that the instability 
of the vortex solution in the attractive case is not driven 
by collapse (^ = eigenvalue) since the vorticity tends to 
push the mass away from the center (r = 0) of the trap. 
This is a general feature of the cases with m > 0, where 
collapse appears to arise in the rings of the cloud, rather 
than at its center. Finally, in Figs. [TTlfT3l we display sim- 
ilar results for the doubly charged case m = 2. 

It is worth mentioning that the unstable modes for 



vorticity-carrying solutions (m > 0) tend to rotate as 
they are growing. This effect can be clearly seen in 
the three-dimensional iso-density contour plot depicted 
in the top panel of Fig. [TH corresponding to the first ex- 
cited state with m = 2 and a = All modes tend 
to stop rotating after the spikes created reach a certain 
height, as can be seen clearly in the bottom panel in 
Fig. [TH This figure shows the ground state with m = 1 
and a = —1. Interestingly, a single solution with multi- 
ple radial nodes (excited states) is able to pick out more 
than one growing mode since each ring can be affected 
by a different g-mode. This effect is seen in the top- left 
panels in Fig. [12] where the first excited state for m = 2 
and cr = —1 is seen to develop, at earlier times, the q = 4 
mode in the inner ring while, at later times, the outer 
ring develops the q = 6 mode. 



IV. CONCLUSIONS 

In this paper, we have revisited the topic of nonlin- 
ear continuation of linear, radially-dependent Laguerre- 
Gauss states of the two-dimensional quantum har- 
monic oscillator model in the presence of inter-particle 
interaction-induced nonlinearity. We have systemati- 
cally constructed such solutions starting from the lin- 
ear limit and, more importantly, we have detailed their 
linear and nonlinear stability properties. This was ac- 
complished by careful examination of the corresponding 
eigenvalue problem or, more appropriately, the (infinite) 
one-parameter family of eigenvalue problems, in radial 
coordinates. We have also provided a full numerical time 
evolution of the model on a radial-polar grid. We have 
principally observed that the ground state of the attrac- 
tive case is unstable due to collapse, although the growth 
rate of the instability may be weak. For the repulsive case 
the relevant state is stable for topological charges m = 
and m = 1; for higher charges the stability depends on 
the atom number. Excited states have been found to 
be generically unstable in both attractive and repulsive 
cases; the development of the instabilities produces col- 
lapse in the former, while it results in the formation of 
azimuthally modulated states in the latter. 

It would certainly be of interest to extend the present 
techniques to the full 3D problem, again considering ra- 
dial states and their continuation from the linear limit, 
as well as their linear stability. However, in the latter 
setting direct numerical computations are significantly 
more intensive. Such studies are outside of the scope 
of the present work and will be considered in a future 
publication. 
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Appendix: Spectral Methods 

When dealing with the Laplacian in polar coordinates, 
the origin presents a significant challenge due to division 
by zero: 



Au 



o u 



1 du 
r dr 



i a u 



(7) 



This problem has been faced with in the past [17|, 118|, \1£ 
[2Q| . with respect to the CP equation, but the methods 
used have so far limited the scope of such investigations. 
However, spectral methods can be used to avoid the r = 
singularity to handle the Laplacian in studies of the CP 
equation. 

The particular application of spectral methods we used 
in our calculations is described by Trefethen in Ref. [13] , 
but is restated here for completeness. This method avoids 
the problem with the origin by avoiding the origin alto- 
gether through the use of an even number of Chebyshev 
nodes. In order to check for the accuracy of utilizing 
spectral methods for the purposes of modeling the CP 
equation in polar coordinates, we recreated the results 
of Fig. lb in Ref. [2^. As can be seen from Fig. [151 
our spectral method code has generated a plot in close 
agreement with the original plot. 

Spectrally discretizing the polar domain is achieved by 
discretizing the radial direction r using polynomial inter- 
polation, in particular the Chebyshev nodes: 



cos(j7r/^) J = 0,1, 



(8) 



meaning the domain must be normalized in order for the 
solution to be properly contained within r = [0, 1], while 
the angular direction 6 is discretized using Fourier (or 
trigonometric) interpolation: 



0j = 2j7r/N j = 0, 1, 



,iV. 



(9) 



These choices for discretization are based upon the pe- 
riodicity of the angular direction and the lack thereof in 
the radial direction (which causes the Gibbs phenomenon 
to occur if Fourier interpolation is used). Additionally, 
the Chebyshev polynomials 



T/e(x) = C0S(/C cos ^{x))^ 



(10) 



can be thought of as a cosine Fourier series {x = cosz), 
meaning they can easily be shown to possess similar re- 
sults for accuracy and convergence as the Fourier method 

[23, m. 

Implementation is simply done by using the nodes de- 
scribed above along with their corresponding derivative 



2000 3000 



FIG. 15: Recreated plot of Figure lb in Ref. [2^ using spec- 
tral methods. The figure shows the unstable part of the main 
eigenvalue as a function of the nonlinear strength N x U for 
a doubly quantized vortex (m = 2) . 



matrices. For the Chebyshev nodes, the Chebyshev dif- 
ferentiation matrix is given by: 



D 



N 



[A, 



Cj {Xi — Xj) 



2{l-x^^ 



2iV2 



6 

2N^ ■ 



1,...,7V-1, 
0, 



6 



where 



Ck 



2 k = 0,N, 
1 otherwise. 



(11) 



(12) 



while the Fourier differentiation matrix (for an even num- 
ber of Fourier nodes) is given by: 



D 




N 



U -i)=0 (modiV) 

cot((j-z)|) (j -i) ^0 (modTV). 

(13) 

It should be noted that the differentiation matrix for an 
odd number of Fourier nodes differs from the one above 
due to a correction which must be made in the deriva- 
tion of the differentiation matrix for an even number of 
Fourier nodes [27], but we have chosen to restrict our- 
selves to using even numbers of Fourier nodes. 

Typically when using polar coordinates, the domain 
of the radial direction is restricted to r G [0, oo] since 
any function on the unit disc must inherently satisfy the 
symmetry condition: 



u{r, 0) = u{-r, {9 + 7r)(mod27r)), 



(14) 
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and since we have chosen to use an even number of 
Fourier nodes, our grid must also satisfy this condition. 
As a result, solution values would be recorded and eval- 
uated twice using this grid. Initially, this may appear to 
be a significant disadvantage for this discretization, but 
an elegant simplification can be arrived at using a sym- 
metry property of the Chebyshev spectral differentiation 
matrix: DN{i, j) = —Dn{N — N — j). 

Using both symmetries, the derivative of u in the radial 
direction can be shown to be exactly the same for both 
occurrences of the node: 

N 
N 

= -DN{N-i, N-j)u{-rj, {0 + 7r)(mod27r)) 

= -<(-r^-, (<9 + 7r)(mod27r)). (15) 

In addition, the derivative calculation can be restricted 
to using only the positive r-axis by using (assuming > 
0): 

N 

K{ri,e) =YDN{iJ)u{rj,0) 
j=o 

N-l 
2 

= DN{iJ)u{rj,e) 

N 

+ Y DN{iJMrj,0) 

N-l 
2 

N 

+ Y DN{iJM-rj, {0 + 7r)(mod27r)) 

N-l 
2 

j=0 

N-l 
2 



The same simplification can be inductively shown to 
work for higher derivatives since they are calculated by 
multiple applications of the derivative matrix 

D2U D'^u D{Du), 

and Eq. (p!5|) has already established the symmetry of the 
first derivative, thus by induction, all higher derivatives 
can be calculated using only the positive r-axis. Conse- 
quently, any simulations based upon a Chebyshev- Fourier 
polar grid can be written to strictly use the positive r- 
axis, resolving the computational time and storage prob- 
lems, but still benefiting from the advantages of spectral 
methods. (See Chapter 11 of Ref. [13] for an implemen- 
tation of this method in Matlab.) 

Spectral methods can also be used to evaluate inte- 
grals with greater accuracy. In particular, this was used 
for calculating the power of the numerically determined 
steady state solutions. To further explore how spectral 
methods can be used for integration, take a generic inte- 
gral 

X 

I{x) = j f{y)dy, (16) 

-1 

for some sufficiently smooth function /. This integral can 
then be restated in the form of an initial value ODE 

^=/(x), /(-1)=0, (17) 

and discretized using Chebyshev nodes, thus replacing 
the derivative by the Chebyshev spectral differentiation 
matrix Dn 

DnI = /, (18) 

where / here is the discrete version of the integrand 
above, / = (/(xq), / (xat-i))"^, and the boundary con- 
dition, /(—I) = 0, is imposed by stripping off the last row 
and column of Dn to obtain the N x N matrix Djy. The 
value of the integral at every Chebyshev node is then 
easily found by inverting the derivative matrix to get 
/ = f . However, the only value we are interested 
in is 1(1) (i.e. the integral over the entire domain) and 
this can easily be found by using only the first row of 
Dn-, call it w\ 1(1) = wf. 
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